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Ubersicht 

Es wird ein neuentwickeltes Finite-Volumen Verfahren fur die dreidimensionalen 
Navier-Stokesschen Gleichungen vorgestellt. Es beruht auf einem Zelleckpunktschema 
mit zentralen Differenzen und expliziten Runge-Kutta Zeitschritten. Durch die Verwen- 
dung lokaler Zeitschritte, impliziter Glattung der Residuen, eines Mehrgitteralgorithmus 
und sorgfaltig kontrollierter kiinstlicher dissipativer Terme wird eine gute Konvergenz 
des Verfahrens zur stationaren Losung erreicht. Es werden Rechenergebnisse fur 
transsonische Profile und Tragflugel angegeben. Es zeigt sich, dap das neue Verfahren 
die roulinemassige Losung der Navier-Stokesschen Gleichungen ermoglicht. 

Einleitung 

Das Gebiet der numerischen Stromungsmechanik hat sich innerhalb der letzten 
Jahrzehnte von der Behandlung einfacher Modellprobleme ausgehend zu einem Stand 
entwickelt, auf dem die numerischen Verfahren bereits eine wichtige Rolle bei der 
Verwirklichung technischer Projekte spielen. Der groPe Fortschritt auf dem Gebiet der 
EDV-Rechenanlagen ermoglichte die Entwicklung einer Vielzahl von numerischen 
Verfahren zur Losung der Eulerschen und Navier-Stokesschen Gleichungen. Die 
Losung der dreidimensionalen Navier-Stokesschen Gleichungen allerdings stellt ein 
Problem dar, fur dessen Bewaltigung die Leistungsfahigkeit der heutigen Vek- 
torrechner gerade ausreicht. Um eine routinema’Pige Anwendung von Losungsverfahren 
fiir diese Gleichungen zu ermoglichen, ist eine Beschleunigung der Verfahren erforder- 
lich. Im vorliegenden Beitrag soli gezeigt werden, dap durch Anwendung der Mehrgit- 
tertechnik in einem Zeitschrittverfahren fiir die Navier-Stokesschen Gleichungen eine 
erhebliche Beschleunigung der Konvergenz zur stationaren Losung erreicht wird. Die 
Leistungsfahigkeit des neuen Rechenverfahrens wird anhand der Rechenergebnisse fiir 
ebene Profilumstromungen sowie fiir Stromungen um den ONERA-M6 Tragfiiigel 
demonstriert. 

Grundgleichungen 

Ausgangspunkt fiir das Verfahren sind die Bewegungsgleichungen fiir die dreidimen- 
sionale, reibungsbehaftete Stromung eines idealen Gases in integraler Form. Sie lauten 
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( 1 ) 


■|-ffJ$dV + J fFftdS = 0 

dt V 6V 

mit als dem Vektor der abhangigen Variablen 

^ = (p pu pv pw pE) T . 


Dabei bezeicnet V ein beliebiges Kontrollvolumen mit der Berandung 5V und der 
auperen Normalen if. F ist der Tensor der FluPdichte nach [1] und die GroPen p, u, v, 
w, E bezeichnen Dichte, Geschwindigkeitskomponenten und Gesamtenergie. Die in F 
enthaltenen Gropen Druck und Temperatur werden mit Hilfe der allgemeinen Zu- 
standsgleichung berechnet. Fur die Berechnung turbulenter Stromungen wird die in den 
Reibungstermen auftretende molekulare Viskositat p durch p+p t und die 
Warmeleitfahigkeit p/Pr durch p/Pr+p^Pr; ersetzt, wobei die Wirbelviskositat p t und 
die turbulente Prandtlzahl Pr t mit einem Turbulenzmodell berechnet werden. In der 
vorliegenden Arbeit wird hierfiir das aigebraische Modell nach Baldwin und Lomax [2] 
verwendet. 


Diskretisierung im Raum 

Die Diskretisierung des Gleichungssystems (1) erfolgt mit Hilfe der Linienmethode, 
d.h. Diskretisierung im Raum und in der Zeit werden getrennt voneinander 
durchgefiihrt. Fur die raumliche Diskretisierung wird das Rechengebiet um den aero- 
dynamischen Koiper in finite Volumen in Form von Hexaedem aufgeteilt, siehe Bild 1. 
Die diskreten Werte der StromungsgrdPen werden den Eckpunkten der Volumen 
zugeordnet. Jeder der Eckpunkte (i,j,k) ist von acht Volumen umgeben. Bild 2 zeigt 
die Anordnung fur den vereinfachten Fall eines ebenen Netzes. Die raumliche Diskre- 
tisierung von Glg.(l) ergibt die diskrete Naherung 



Im folgenden sollen die diskreten reibungslosen Fliisse (3c jjk. die reibungsbehafteten 
Fliisse (3 d y t k sowie die kiinstlichen dissipativen Terme erlautert werden. 

Fur die Bilanz der reibungslosen Fliisse wird die Zelle ABCDEFGH nach Bild 2 
betrachtet. dc ijk durch Summation der Fliisse fiber die Berandungsflachen der 
Zelle ermittelt. Die FluPdichte auf den Randflachen wird durch Mittelwenbildung der 
GroPen auf den benachbarten Knoten gebildet und mit dem diskreten Flachenvektor 
multipliziert. Es la'Pt sich zeigen [3], dap sich der Diskretisierungsfehler wie von erster 
Ordnung verhalt, wenn der Verlauf des Normalen vektors fiber den Randflachen der 
Zelle glatt ist und die Form der Randflachen sich bei Netzverfeinerung einem Paral- 
lelogramm nahert. Auf glatten Rechennetzen verhalt sich der Diskretisierungsfehler 
wie von zweiter Ordnung. 

Ffir die Bilanz der reibungsbehafteten Fliisse wird die Hilfszelle A B C D betrachtet. 
Die Fliisse do i.j.k enthalten erste Ableitungen der Stromungsgropen und werden im 
Mittelpunkt der Berandungsflache durch Anwendung einer lokalen Transformation und 
finite Differenzen approximiert. Als Alternative bote sich die Ermittlung der ersten 



Ableitungen in mit Hilfe des Greenschen Satzes und die Berechnung der 

FluPgr6(3en durch Mittelung an. Diese Vorgehensweise fiihrt nach [4] zu einer Entkop- 
plung benachbaner Gitterpunkte Oder sogar zur Anfachung von Oszillationen in der 
Losung je nach dem Seitenverhaltnis des Hexaeders. Fiir den Fall eines eindimen- 
sionalen gestreckten Netzes la(3t sich zeigen [4,5], daP sich der Diskretisicrungsfehler 
der reibungsbehafteten FluPbilanz wie von erster Ordnung auf beliebig gestreckten 
Netzen und wie von zweiter Ordnung auf glatten gestreckten Netzen verhalt. Im vor- 
liegenden Verfahren wurden bei der Ermittlung der reibungsbehafteten Fliisse nur Gra- 
dienten in der Koordinatenrichtung tj normal zum Korper beriicksichtigt (thin layer 
approximation). 

Die kiinstlichen dissipativen Terme wirken dampfend auf Oszillationen im 

Stromungsfeld und ermoglichen somit die Konvergenz des Verfahrens zur stationaren 
Losung. Sie bestehen aus einer Kombination zweiter und vierter Differenzen. Um fur 
Volumen mit hohen Seitenverhaltnissen, die bei der Berechnung von reibungsbehaf- 
teten Stromungen zwangslaufig auftreten, gute Dampfungseigenschaften zu erreichen, 
muP die Grbpe dieser Terme sorgfaltig festgelegt werden. Wir verwenden 

% = (D\ - T>\ + - D\ + D\ - D\) lik (3) 

mit den zweiten und vierten Differenzenoperatoren 

W ij,k = i+l/2,j,k' e(2) i+l/2.j.k)^^i.j.k > ( 4 ) 

D 4 5 vfio.t = , (5) 

wobei ij,k die Indizes in £,T|,£-Richtung bedeuten und A^,V^ vorwarts- und 
riickwartsgerichtete Differenzenoperatoren in ^-Richtung sind. In Anlehnung an Mar- 
tinelli [6] ist der Skalierungsfaktor vom grbpten Eigenwert der Jakobimatrix der 
reibungslosen Fliisse 3(^ c /3$ sowie dem Seitenverhaltnis des Volumens abhangig: 

\ i+l/2j,k “ \ i+l/2j,k ‘ ^i+ltfj.k (6) 

^i+l^.j.k = 1 + max ((Xq i+l/2.j,k) a » (X; i+l/2j,k / X£ | i+l/2,j,k) a ) 

X' i+l/2,j,k = l^i+l/2,j,k * ^ i+l/2,j,k * + C 1^ i + ]/2,j,k 1 

Xq i+l/2,j,k = ^i+l/2,j,k * ^>q i+l/2,j,k I + c ^q i+l/2.j,k I 

X; i+l/2j.k = '^i+]/2.j.k ’ ^q i+l/2j.k I + C 1?^ j + i/2.j,k 1 

wobei tfi+i/2,j,k ber Geschwindigkeitsvektor und j+j^j.k* ^q i+i/2.j,k’ i+i/2,j,k die 

diskreten Zellflachenvektoren in J;,r|,£-Richtung und c die Schallgeschwindigkeit 
bedeuten. Der Exponent a ist l/2<a<2/3. Die Koeffizienten e^und e (4 ^ verwenden den 
Druck als Sensor fiir VerdichtungsstoPe im Stromungsfeld, sie lauten 

e(2) i+l/2,j,k = k (2) max( v i_i j k . v i,j.k. v i+l j,k» v i+2,j,k) (X) 

_ | Pi-l.j.k ~ Xpj j.k + Pj+l.j.k | 

Pi-l,j,k XpjJ.k + Pi+],j,k 

e (4) i + i/2j,k = max ( 0 • ( k(4) " e (2) i + i/2.j.k)) 





mit den Konstanten k® und k^. Die Dissipationsoperatoren in T]- und ^-Richtung sind 
cntsprechend definiert. 

Randbedingungen 

An den Femfeldrandem des Rechengebietes werden die Strbmungsgrb|3en entsprechend 
dem Verlauf der Charakteristiken fur eindimensionale Stomung normal zum Rand nach 
[7] ermittelt. Auf der Korperoberflache wird die Geschwindigkeit zu null gesetzt. Die 
Kontinuitatsgleichung und die Energiegleichung werden fur die Gitterpunkte auf der 
Korperkontur unter Annahme einer adiabaten Wand gelost. 

Zeitschrittverfahren , 

Durch die Diskretisierung im Raum erhalt man ein System gewohnlicher 
Differentialgleichungen, die mit einem funfstufigen Runge-Kutta Verfahren gelost wer- 
den, Entsprechend [6] verwenden wir ein Schema, bei dem die kiinstlichen dissipativen 
Terme in den Stufen eins, drei und fiinf berechnet werden. Dieses Schema hat einen 
besonders weiten Stabilit'atsbereich beziiglich der dissipativen Terme, was bei Rechen- 
netzen mit stark variierendem Seitenverhaltnis der Zellen von Bedeutung ist. Die 
stationare Losung des Gleichungssystems ist unabhangig von der Gro(3e des 
Zeitschritts und somit auch von den im folgenden beschriebenen Beschleunigungstech- 
niken (1) der Iokalen Zeitschritte, (2) der impliziten Glattung der Residuen und (3) des 
Mehrgitteralgorithmus. 

Bei der Methode der Iokalen Zeitschritte wird in jedem Gitterpunkt der grofite von der 
Stabilitatsbedingung des expliziten Verfahrens her erlaubte Zeitschritt durchgefiihrt. In 
der Berechnung des Zeitschritts sind die Stabilitatsgrenzen beziiglich Konvektion und 
Diffusion enthalten [1], 

Die expliziten Residuen R^ des Zeitschrittverfahrens werden in jeder Runge-Kutta 
Stufe einer impliziten Glattung unterworfen. Fur dreidimensionale Probleme verwenden 
wir die faktorisierte Form 

(l-£ 5 V t A 4 )(l-^V^)(l-E 5 V ? A 5 ) Rn, = Rn, (8) 

wobei Rj,, das zu berechnende implizite Residuum der m-ten Stufe bedeutet. Die 
Glattung der Residuen ist aus zwei Grunden vorteilhaft. Zum einen erhalt man einen 
vergroperten Stabilitatsbereich des Verfahrens, die Zeitschritte konnen im Vergleich 
zum expliziten Schema etwa verdoppelt werden. Zum anderen kann aber auch das 
Dampfungsverhalten des Verfahrens in Bereichen mit Zellen hohen Seitenverhaltnisses 
verbessert werden. Das Dampfungsverhalten des Runge-Kutta Zeitschrittverfahrens ist 
namlich besonders gunstig, wenn der Zeitschritt nur etwas ldeiner als von der explizi- 
ten Zeitschrittbegrenzung her erlaubt ist. Fur sehr viel kleinere Zeitschritte ver- 
schlechtert sich das Dampfungsverhalten. Zur weiteren Erlauterung nehmen wir an, es 
sei eine hochfrequente Storung in der transienten Stromungslbsung entlang der langen 
Seite einer Rechenzelle vorhanden. Die explizite Zeitschrittbegrenzung wird durch die 
Lange der kurzen Seite der Zelle bestimmt, so dap sich das Verfahren beziiglich des 
Dampfungsverhaltens entlang der langen Seite unterhalb des Optimalpunktes bester 



Dampfung befindet. Durch eine Variation dcr Glattungskoeffizienten e^e^.e^ in 
Abhangigkeit vom Zeitschritt und vom Seitenverhaltnis der Zellen, z. Bsp. durch 


= max(0 , — ( 


CFL 


4 CFL explizil X^+max^,^) 


(l+max((-^-) a > (^-) a ))) 2 -l) 


(9) 


kann das Dampfungsverhalten wesentlich verbessert werden. 


Der Mehrgitteralgorithmus basiert auf dem Prinzip, die Losung auf einem 
feinen Netz unter Zuhilfenahme groberer Rechennetze zum stationaren Zustand zu kon- 
vergieren. Da die groberen Rechennetze weniger Rechenaufwand erfordem und zudem 
gropere Zeitschritte zulassen, kann man eine Beschleunigung des Verfahrens erwarten. 
Femer wird mit Hilfe der groben Netze das Dampfungsverhalten beziiglich der 
langwelligen Storungen im Stromungsfeld verbessert, so dap auf die sonst iibliche 
Enthalpiedampfung des Zeitschrittverfahrens verzichtet werden kann. Der Mehrgit- 
teralgorithmus ist entsprechend der Arbeit von Jameson [8] aufgebaut. Ausgehend vom 
feinsten Rechennetz werden grobere Netze durch Weglassen jedes zweiten Gitter- 
punktes in jeder Koordinatenrichtung erzeugt. Die Stromungsvariablen auf den 
Netzpunkten der groben Netze werden vom jeweils feineren Netz ubemommen. Die 
Residuen der groben Netze werden mit den Residuen der jeweils feineren Netze so 
verkniipft, dap die Losung auf den groben Netzen von den Residuen der feinen Netze 
angetrieben wird, und somit die auf den groben Netzen berechneten Korrekturen der 
Losung null sind, wenn die Residuen des feinsten Netzes null sind. Nach wiederholter 
Vergroberung des Netzes und Erreichen des grobsten Netzes werden dessen Korrek- 
turen zum nachst feineren Netz trilinear interpoliert und zu den Korrekturen dieses 
Netzes addiert. Die Interpolation der Korrekturen zum nachst feineren Netz kann nun 
entweder solange wiederholt werden, bis das feinste Netz erreicht wird (V- 
Mehrgitterzyklus), oder es konnen zwischendurch noch Schritte zuriick auf grobere 
Netze durchgefiihrt werden. Es zeigt sich, dap mit einem W-Mehrgitterzyklus eine 
schnellere Konvergenz zur stationaren Losung erzielt werden kann als mit dem V- 
Zyklus. Fur das Dampfungsverhalten des Verfahrens ist es vorteilhaft, wenn die 
Summe der Korrekturen der groben Netze einer Glattung nach Glg.(8) mit 
6^,6^, = 0.2 unterworfen werden, bevor sie auf die Losung des feinsten Netzes ad- 
diert wird. Die Robustheit des Verfahrens wird ferner verbessert durch sukzessive 
Netzverfeinerung. Hierbei wird das Verfahren zunachst auf einem groben Netz gestar- 
tet und dessen Losung dient als Anfangsbedingung fur das nachst feinere Netz. 


Ergebnisse 

Einige Eigenschaften des neuen Rechenverfahrens sollen zunachst anhand der Ergeb- 
nisse einer vereinfachten Version fur ebene Stromungen dargestellt werden. Als Test- 
fall wahlen wir die transsonische Umstromung des RAE 2822 Profils, 
Ma^^.73, a=2.79°, Re=6.5T0 6 . Es wird ein C-formiges Rechennetz mit 385x65 
Netzpunkten verwendet. Der Abstand des ersten Gitterpunktes von der Wand betragt 
10 -5 Profillangen. Bild 3 zeigt den EinfluP des Dissipationskoeffizienten k (4) auf das 
Konvergenzverhalten des Verfahrens. Die Spitzen im Residuenverlauf zu Beginn der 
Rechnungen zeigen die sukzessive Netzverfeinerung an. Mit einem mittleren Wert von 



k w werden etwa 500-600 Mehrgitterzyklen auf dcm feinsten Netz benorigt, um 
I Idp/dtl l 2 um 10 Gro(3enordnungen zu vennindem. Femer zeigt sich, da(3 der 
Dissipationskoeffizient um nahczu eine GroPenordnung variiert werden kann, ohne dap 
die Konvergenz des Verhahrens verloren geht. Das Verhalten des Verfahrens 
beziiglich der Netzfeinheit und den Vergleich der Rechenergebnisse mit Messungen [9] 
zeigt Bild 4. Die globalen Eigenschaften der Stromung werden bereits mit dem groben 
193x33 Netz aufgelost. Die Unterschiede zwischen den Netzen 385x65 und 577x97 
sind sehr gering. Auffallig ist femer, dap der Verlauf der Wandschubspannung schon 
mit 33 Netzpunkten in Normalenrichtung recht genau berechnet wird. Tabelle 1 zeigt 
anhand verschiedener Konvergenzkriterien, daP eine fur praktische Anwendungen hin- 
reichend konvergierte Losung innerhalb von 100 Mehrgitterzyklen erreicht wird, das 
entspricht bei einem 385x65 Netz einer Rechenzeit von etwa 150 s auf Cray-2. Ohne 
Mehrgitteralgorithmus werden dagegen je nach Konvergenzkriterium zwischen fiinf- 
und zehnfach hohere Rechenzeiten benotigt. 

Als nachstes wird das Rechenverfahren auf die transsonische Umstromung des 
ONERA-M6 Tragfliigels angewendet. Das Rechennetz weist eine C-Struktur im 
Profilschnitt und eine O-Struktur im Spannweitenschnitt auf, Es werden 289 
Netzpunkte im Profilschnitt, 65 normal zur Oberflache und 49 in Spannweitenrichung 
verwendet. Der Abstand des ersten Netzpunktes von der Oberflache betragt 10 -5 der 
ortlichen Fliigeltiefe. Fiir die Anstrombedingungen Ma^^.84, a=3.06°, Re=l 1-10 6 
liegt eine weitgehend anliegende Stromung vor. Bild 5 zeigt den Konvergenzverlauf 
der Rechnung und die Druckverteilungen in verschiedenen Spannweitenschnitten, 
wobei auch die Ergebnisse auf einem etwas groberen Netz mit 193x49x33 Punkten mit 
angegeben sind. Insgesamt liegt eine gute Ubereinstimmung mit den Messungen nach 
[10] vor. Unterschiede in den Ergebnissen Fur die beiden Netze ergeben sich in der 
Nahe von Verdichtungsstopen, wo das feinere Netz die groPen Gradienten der 
StromungsgroPen besser auflost. Fiir • die Anstrombedingungen 
Ma^^.84, a=6.06°, Re=ll-10 6 zeigt Bild 6, dap die Ubereinstimmung zwischen Mes- 
sung und Rechnung hier weniger gut ist. Anhand der berechneten Wandstromlinien 
erkennt man die Bereiche abgeloster Stromung auf der Fliigeloberseite. Das Tur- 
bulenzmodell [2] liefert hier zu hohe Werte der Wirbelviskositat p t , was zu einer zu 
kleinen Abloseblase fiihrt. Aus Tabelle 1 kann entnommen werden, daP Fiir die prak- 
tische Anwendung des Verfahrens 100 Mehrgitterzyklen fiir eine stationare Losung 
ausreichen. Das entspricht bei dem 193x49x33 Netz einer Rechenzeit von 40 min auf 
Cray-2. 

Zusammenfassung 

Es wurde liber ein Finite-Volumen Verfahren fiir die Navier-Stokesschen Gleichungen 
berichtet. Es beruht auf einem Zelleckpunktschema mit zentralen Differenzen und 
einem expliziten mehrstufigen Runge-Kutta Zeitschrittverfahren. Fiir eine gute Kon- 
vergenz zur stationaren Losung werden lokale Zeitschritte, eine implizite Glattung der 
Residuen, ein Mehrgitteralgorithmus und sorgfaltig kontrollierte kiinstliche dissipative 
Terme verwendet. Rechenergebnisse zeigen, daP konvergierte Losungen fiir ebene 



transsonische Profilumstromungen mit 25000 Netzpunkten 2-3 Minuten und fur dreidi- 
mensionale Tragfliigel mit 300000 Punkten 40 Minuten Rechenzeit auf einem Cray-2 
Rechner benotigen. An dem Einbau altemativer Turbulenzmodelle fiir eine genauere 
Berechnung von Stromungen mit gro(3en Ablosegebieten wird derzeit gearbeitet. 

i 
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FaD 

Netz 

k (4) 

1%C A 

0.1% C A 

10' 5 Red. 

a 

■ 

RAE 2822 

385x65 

1/64 

23 (50) 

49 (82) 

278 (346) 

0.8402 

0.01738 

Ma_=0.73 

385x65 

1/32 


49 (82) 

198 (265) 

0.8356 

0.01730 

cc=0.73 

Re=6.5-10 6 

385x65 *) 

1/32 

1080 (518) 

1970 (957) 

2950 (1418) 

0.8356 


0NERA-M6 

Ma„=0.84 

0=3.06° 

Re=11.010 6 

288x65x49 

1/32 

6 (886) 

54 (3745) 

-460(27580) 

0.2689 

0.01784 

0NERA-M6 

Ma*,=0.84 

0=6.06° 

Re=11.010 6 

288x65x49 

1/32 ■ 

21 (2012) 

44 (3362) 


0.5380 

0.05559 


*) ohne Mehrgitteralgorilhmus, nur feinstes Netz 


Tabelle 1 : Konvergenzverhalten des Verfahrens fur zwei- und dreidimensionale 
Stromungen. Zahlen geben die fur das jeweilige Konvergenzkriterium erford erlichen 
Mehrgitterzyklen auf dem feinsten Netz an. Zahlen in Klammem geben die Sekunden 
in Rechenzeit auf Cray-2 an. 



Bild 1 : Dreidimensionales korperangepafhes Rechennetz 
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Bild 3 : EinfluP des Dissipationskoeffizienten k^ auf das Konvergenzverhalten fur 
RAE 2822 Profil, Ma^O.73, a=2.79°, Re=6.510 6 , Netz 385x65 . 
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a) 193x33 Netz b) 384x65 Netz c) 577x97 Netz 

Bild 4 : Druckverteilungen und Wandschubspannung fur verschiedene Netzdichten, 
RAE 2822 Profil, Ma«=0.73, a=2.79°, Re=6.5-10 6 , k (4) =l/64 . 
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log I Idp/Otl I 



Bild 5 : Konvergenzverhalten und Druckverteilungen fur ONERA-M6 Fliigel, 
Ma^.84, oc=3.06°, Re=ll-10 6 . 
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log I I dp/3t I I 



Bild 6 : Konvergenzverhalten, Druckverteilungen und Wandstromlinien fiir ONERA- 
M6 Flugel, Ma.^.84, oc=6.06°, Re=lM0 6 . 




